Für die meisten GIS Anwendungen in R könnt ihr das package “sf” benutezen. Eine ausführliche Anleitung zu fast allem, was sf zu bieten hat findet ihr hier. Im Rahmen des Projektes werdet ihr sf brauchen, um die einzelnen Probelstellen ihren Catchments zuzuweisen. In diesem Skript werden wir die Probestellen Bundesländern zuordnen, und uns die Catchments auf Detuschland beschränken.
Als erses laden wir das Package.
library(sf)
## Linking to GEOS 3.9.1, GDAL 3.2.1, PROJ 7.2.1
Dann zuerst die Probendaten. Die Probedaten liegen im .RDS Format vor. Dieses Format ist speziell für R und alle R Objkete lassen schnell darin speichern saveRDS() und landen readRDS(). Besonders bei rumlichen Daten kann die größe der Dateien und die Ladedauer damit stark reduziert werden. Die Daten können allerdings von keinem anderen Programm (z.B. QGIS) geöffnet werden.
mzb <- readRDS("data/germany.rds")
und die Catchments. Die Catchments liegen als Shapefile vor. Shapefiles kennt ihr hoffentlich gut aus dem GIS Kurs. Sie bestehen immer aus mehreren Datein aber nur die .shp Datei kann genutzt werden um die räumlichen Daten zu öffnen. Um die Shapefile in R einzulesen nutzen wir die erste Funktion aus sf: sf_read. Alle Funktionen in sf fangen mit dem beiden Buchstaben st (kurz für Spatio-temporal) an. Da R Studio über eine autocomplete Funktion verfügt kann man also immer sf_ eintippen und bekommt die vollständige Liste an Funktionen aus dem sf Package.
cat <- st_read("data/lemm/MultipleStress_RiverEcoStatus.shp")
## Reading layer `MultipleStress_RiverEcoStatus' from data source
## `C:\Users\jonat\Documents\Uni\teaching\projekt_uwi\21_betadiversity\data\lemm\MultipleStress_RiverEcoStatus.shp'
## using driver `ESRI Shapefile'
## Simple feature collection with 52847 features and 11 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: -10.4779 ymin: 34.64526 xmax: 32.93388 ymax: 70.1431
## Geodetic CRS: WGS 84
Einen Polygonlayer der deutschen Bundesländer können wir mit den geodata package runterladen. Den Code werde ich hier ich explizit erklären, mehr Infos zu dem geodate package findet ihr hier.
library(geodata)
## Lade nötiges Paket: terra
## terra version 1.4.20
ger1 <- gadm(country = "DEU", level = 1, path = "data")
ger1 <- st_as_sf(ger1)
Im folgenden Abschnitt werde ich einge Funktionen demonstrieren die man allgemein häufig gebrauchen kann. Die ersten Funktionen (filter(), select(), ) stammen aus dem dplyr package. Man benutzt sie um bestimmte Reihen (filter) oder Spalten (select) einer Tabelle auszuwählen. Die Probedaten habem viele Spalten (n = ncol(mzb)), die meisten davon sind für uns hier uninterissant. Wir wählen nur die site_id, den Genusnamen und das Datum aus. Wenn man select auf ein sf Objekt anwendet wird die Spalten mit den Koordinaten automatisch mit ausgewählt.
library(dplyr)
##
## Attache Paket: 'dplyr'
## Die folgenden Objekte sind maskiert von 'package:terra':
##
## intersect, src, union
## Die folgenden Objekte sind maskiert von 'package:stats':
##
## filter, lag
## Die folgenden Objekte sind maskiert von 'package:base':
##
## intersect, setdiff, setequal, union
mzb2 <- select(mzb, c("site_id", "genus", "date"))
head(mzb2)
## site_id genus date
## 1 00001 Antocha 2014-04-15
## 2 00001 Ancylus 2014-04-15
## 3 00001 Baetis 2014-04-15
## 4 00001 <NA> 2014-04-15
## 5 00001 Dicranota 2014-04-15
## 6 00001 Dugesia 2014-04-15
Aus diesem Datensatz können wir jetzt noch die Reihen auswählen die Baetis Arten enthalten. Wie immer in R müssen wir zwei Gleichzeichen benutzen um ist gleich auszudrücken, da ein einfaches Gleichzeichen auch für Objektzuweisungen genutzt werden kann.
mzb3 <- mzb2 |> filter(genus == "Baetis")
Anbei ein weiteres beispiel in dem wir drei Bundsländer aus dem ger1 Objekt herrausfiltern. Das %in% Zeichen bedeuted ist Element von. Der Befehl sucht also alle Reihen aus ger1 aus in dennen die Spalte NAME_1 einem der Elemente des Vektors `c(“Bayern”, “Berlin”, “Bremen”)``entspricht.
ger2 <- filter(ger1, NAME_1 %in% c("Bayern", "Berlin", "Bremen"))
Das können wir uns auch auf einer Karte ansehen:
library(mapview)
mapview(ger2)
Manchmal möchte man auch nur eine Spalte entfernen oder alle Reihen auswählen die einen bestimmten Eintrag nicht haben. Das geht in R mit dem Ausrufezeichen. Beispiel:
cat2 <- select(cat, !eco_stat_n)
cat2 <- filter(cat, eco_stat_2 != "bad")
cat2 <- filter(cat, !eco_stat_2 %in% c("bad", "moderate"))
Die Probedaten sind aktuell kein räumliches Objekt. Sie sind eine einfacher data.frame (eventuell wird bei euch data.table angezeigt). Um aber Probedaten auszuwählen, die in einem bestimmten Bundesland oder Catchment liegen müssen sie räumlich werden - also sf Objekte. Das machen wir mit der Funktion st_as_sf(). Dabei gibt es drei Sachen zu beachten. 1. Aktuell entspricht jede Reihe in mzb der Beobachtung von einer Art, an einem Datum und an einer Probestelle. Würden wir diese Tabelle so in eine sf Objekt umwandeln, hätten wir für jede Art and jeden Datum einen Punkt. Da viele der Punkte übereinander liegen würden bringt uns das nichts. Wir können den Rechenaufwand also reduzieren, indem wir nur eine Reihe pro Probestelle behalten. Auch dafür können wir eine Funktion aus dplyr nutzen. Mit distinct() wählen wir eine Reihe pro Eintrag der angegebenen Variable (site_id) aus.
mzb2 <- distinct_(mzb, "site_id", .keep_all = TRUE)
## Warning: `distinct_()` was deprecated in dplyr 0.7.0.
## Please use `distinct()` instead.
## See vignette('programming') for more help
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was generated.
mzb$EPSG aus da die Funktion verwirtt ist wenn wir einen Vektor als Argument für das Koordinatenbezugssystem angeben.sf_mzb <- st_as_sf(mzb2, coords = c("x.coord", "y.coord"), crs = mzb$EPSG[1])
Letztlich wollen wir hier, als Beispiel, nur die Probestellen raussuchen, die in Bayern liegen. Wir müssen darauf achten, dass beide Objekte (die Probestellen und Bayern) im gleichen KBS liegen. Das können wir mit der Funktion st_crs prüfen.
st_crs(sf_mzb) == st_crs(ger1)
## [1] FALSE
Aktuell haben sie noch unterschiedliche KBS. Wir müssen also einen von beiden transformieren.
ger2 <- st_transform(ger1, crs = st_crs(sf_mzb))
Die Auswahl der Probestellen erfolgt durch st_filter(). Die Funktion filtert alle Objekte aus dem ersten Argument die innerhalb des zweiten Argumentes liegen. Ganz ähnlich funktioniert die Funktion st_join() mit der ihr Variablen von einem Layer (den Catchments) auf einen anderen (die Probestellen) übertragen könnt.
bav <- filter(ger2, NAME_1 == "Bayern" )
sf_mzb_bav <- st_filter(sf_mzb, bav)
mapview(sf_mzb_bav)
Für die Analysen würde es Sinn machen sich auf die aktivsten Jahre (1986:2013) und die Monate April bis Juni zu konzentrieren.
Um aus einem Datum, wie der date Spalte im Probedatensatz den Monat zu extrahieren könnt ihr die Funktion month() aus dem lubridate package verwenden.
library(lubridate)
##
## Attache Paket: 'lubridate'
## Die folgenden Objekte sind maskiert von 'package:terra':
##
## intersect, union
## Die folgenden Objekte sind maskiert von 'package:base':
##
## date, intersect, setdiff, union
months <- month(sf_mzb_bav$date)
Die Polygone der Shapefiles sind leider nicht gut gemacht. Was genau das Problem ist weis ich nicht, es könnte sein, das einige Polygone nicht perfekt schließen oder sich die Linien eines Polygons selber kreuzen. Wie dem auch sein all das können wir mit der wunderbaren Funktion st_make_valid() beheben.
cat2 <- st_make_valid(cat)
Wie oben erwähnt wollen wir nur die einzelnen Catchments die aktuell noch große Teile Europas bedecken auf Detuschland reduzieren. Das geht mit der bereits benutzen st_filter() Function.
cat2 <- st_transform(cat2, crs = st_crs(ger1))
cat3 <- st_filter(cat2, ger1)
mapview(cat3)